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SECTION  1 


INTRODUCTION 

1.1  INTRODUCTION 

The  purpose  of  this  research  is  to  investigate  diffraction  effects  in 
complicated  geologic  and  topographic  structures  of  interest  to  the  Air  Force 
Geophysics  Laboratory.  This  first  technical  report  examines  diffraction  by 
reentrant  corners,  i.e.,  corners  with  included  angle  greater  than  180  degrees, 
which  idealize  the  abrupt  topographic  relie*  commonly  found  at  the  base  of 
escarpments  and  mesas.  Reentrant  landforms  are  also  common  in  tectonic 
regions,  and  include  fault  scarps  and  other  block  joints.  These  erosiona! 
and  tectonic  features  are  nicely  illustrated  in  Fig.  1,  reproduced  collectively 
from  Holmes  (1965),  showing  near-surface  structural  variations  at  both  local 
and  regional  scales.  Approximating  these  regions  of  high  curvature  by  perfect 
edges  defines  so-called  generalized  model s--where  curvature  is  concentrated 
at  singular  points  (corners)  on  the  free-surface.  Diffraction  from  a  single 
corner  will  be  studied  in  detail  here.  Effects  of  more  realistic  distributed 
curvature  in  local  and  regional  models  will  be  examined  in  subsequent  reports. 

The  principal  tool  for  the  study  is  finite  element  modeling,  since  no 
mathematical  solutions  presently  exist  for  P-SV  wave  edge  d4 'fraction  in 
even  the  simplest  geologic  d' ffractors .  This  necessary  emphasis  on  numerical 
simulation  is  balanced  in  the  study  by  a  mathematical  analysis  of  self¬ 
similar  vector  wave  diffraction  theory.  The  analysis  yields  a  formal  reduc¬ 
tion  of  the  general  problem  to  a  system  of  singular  integral  equations.  The 
resulting  integral  expressions  provide  a  new  calculational  basis  for  solutions 
to  a  variety  of  diffraction  orob7ems  in  mathematical  physics. 

1.2  BACKGROUND 

Motivation  *or  this  study  is  the  observation  that  highland  interiors 
are  often  in  the  shadow  zone  o*  low’and  seismic  sources.  Consequently , 
diffracted  waves  surfaces  of  high  curvature  (e.g.,  edges  and  corners) 
provide  the  dominant  seismic  s^Gne’s  detected  in  the  interior,  in  other 
words,  edges  and  corners  on  tl,e  *ree  sur*ace  act  as  secondary  sources 
i 1  ’ u"i nation  tk-o  VeHland  interior.  This  is  4 1  lustrated  in  Fig.  2,  showing 
a  cartoon  cr  tHe  gnlard- lowland  tcpograohy,  source  generated  rays  and  the 
geometric  shadow  rc-’e.  arc  d*’f  fraction  (secondary)  sources.  Numerical  cal  - 
cu’at'on-  r'r  t"'1  * wave  r4°’  d  are  "ecessarv  to  simulate  this  phenomenon 


because  approximate  theories,  such  as  conventional  ray  tracing  and  boundary 
integral  methods,  cannot  accommodate  diffraction  from  irregular  boundaries. 
However,  questions  remain  regarding  the  range  of  validity  of  numerical  dif¬ 
fractions,  due  in  large  part  to  the  existence  of  certain  singularities  (stress 
concentrations)  at  reentrant  corners  (geometric  singularities). 

The  geophysical  problem  is  best  illustrated  by  some  finite  element  ex¬ 
amples  modeling  reentrant  wedges  of  elastic  rock  diffracting  an  incident 
Rayleigh  wavelet.  Snapshots  of  the  full  wave  field  are  shown  in  Fig.  3,  for 
Rayleigh  wave  incidence  from  the  right.  Each  wedge  simulates  the  lower 
edge  of  an  escarpment,  and  the  Rayleigh  wave  represents  ground  roll  from  a 
vibratory  surface  source  located  100  meters  (100  elements)  from  the  edge. 
Although  the  escarpment  is  in  the  source's  geometric  shadow,  clearly  surface 
wave  interaction  with  the  corner  radiates  significant  body  wave  energy  into 
the  highland  interior.  Furthermore,  the  Rayleigh  wavelet  is  seen  to  transmit 
effectively  around  the  corner,  hence,  will  radiate  into  the  interior  at  each 
encounter  with  subsequent  corners  on  the  free  surface. 

1.3  REPORT  PLAN 

Results  of  the  research  are  presented  in  the  following  three  sections. 
Section  2  describes  the  numerical  diffraction  experiments,  Section  3  de¬ 
velops  the  mathematical  theory  of  vector  wave  diffraction  and  Section  4 
presents  the  discussion  and  conclusions.  The  numerical  work  in  Section  2  was 
done  principally  by  the  second  author.  The  first  author  was  responsible  for 
the  remaining  research  and  written  report. 


Fig.  791  Schematic  section  across  the  Cordillera  of  the  western  United  States 
Length  of  section  about  !,900  miles  (3,000  km.) 


Figure  1-1.  Examoles  of  erosionai  and  tectonic  Tandems  with  reentrant 

featu>-es.  The  upper  two  diagrams  are  at  local  scales  (l-100km) 
while  the  bottom  section  is  on  a  regional  scale  (1000-3000km) . 
~hese  illustrations  with  cautions  are  reproduced  directly  from 
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1-2.  Illustrations  of  a)  highland-lowland  erosional 
topography  showing  a  typical  mesa  cross-section; 
b)  idealized  topography  showing  rays  from  a  sur¬ 
face  source  on  the  lowland  and  its  highland  geo¬ 
metric  shadow  zone;  c)  diffractions  from  the  cor 
ners  illuminating  receivers  in  the  shadow  zone. 


sequence.-  or'  snapshots  from  finite  element  simu- 
o  tiers  of  Rayleigh  Surface  wave  reflection,  trans 
ission  ana  diffraction  m  four  reentrant  corners, 
he  Rayleigh  wave  is  incident  from  the  right,  and 
ach  sequence  shows  vector  plots  of  the  velocity 
ie'. d  at  50  times tep  increments,  starting  when  the 
jrface  wave  reaches  the  corner  at  timestep  =  250 


r.u...er  ta .  _.f  c  .fractio.-i  Trom  reentrant  corners  is  performed 

uS.r.g  i.'.c  ri.iite  e.er.erit  7.c...Cv..  Aiane  A- ,  oV-  anc  Pri-wavei  ano  —  Aayieign  sur- 
"iCi  Ac..i  _re  i r.ci cent  on  reentrant  corners  with  angles  typical  of  mesas  ana 
escarpments.  Vr.e  Sn-waves  have  partic'. e  motions  perpendicular  to  the  plane  of 
cr.c  met  a  itr^ct-re,  1  e  ?-SV  Wcves  nave  ir.-p'.ane  particle  motions.  ~r,e  Ray¬ 
leigh  waves  r.e^  c  in- plane  particle  motions  witn  their  cnaracteristic  retrograde 
elliptical  oroits  ana  exponential  decay  with  depth.  The  resulting  diffracted 
fields  are  then  examined  quantitatively. 

Section  2.2  outlines  the  finite  element  model  used  for  the  computations. 
Section  2.3  discusses  diffractions  for  SH-wave  input.  Results  for  P-SV-wave 
'•.‘.put  are  snowt  in  Section  2.4, while  Section  2.5  concludes  with  the  results 
for  Rayleigh  wave  input. 

2.2  FINITE  ELEMENT  MODEL 

Tr.e  numerical  computations  are  made  with  FLEX  ,  an  explicit  time  in¬ 
tegration  finite  element  code,  Vaughan  (1984).  Tne  baseline  finite  element 
•'•'-del  is  a  i99  meter  cv  199  meter  cross-section  of  a  mesa  edge  with  a  199  by 
.99  element  ci bcretization  of  1  meter  square  elements,  Fig.  2-1.  The  material 
ii  linearly  elastic  with  c^  =  2600  m/sec,  c$  =  1500  m/sec  and  p  =  2000  kg/m3. 
Frequency  resolution  of  the  model  is  aoout  160  Hz.  Reentrant  angles  are 
mcaeiea  by  cutouts  in  the  first  and  second  quadrants  with  the  vertex  at  the 
center  of  tne  gric.  Note  that  definition  of  the  mesa  by  square  elements  pre¬ 
cedes  a  steppeo  interface  for  most  of  the  reentrant  angles  considered. 

Body  waves  are  input  from  the  bottom  of  the  grid,  the  left  and  right  siae^ 
.saving  appropriate  symmetry  or  asymmetry  boundary  conditions  for  P-  or  Sv-waves 
respectively.  Rayleign  surface  waves  are  generated  by  prescribing  vertical 
velocity  over  five  nodes  (four  elements)  along  the  x-axis.  Fig.  2-i,  with 
a  symmetry  condition  on  the  rignt  side  and  an  absorbing  boundary  on  the 
other  three  sides.  All  wave  types  are  input  as  nodal  velocity-time  his¬ 
tones  of  tr.e  wavelet  function,  f  -  A*EXP(-  (uT/a)2)  cos(JT  +  v ) ,  where  A  is 
amp. itude, ~ is  predominant  wavelet  frequency , a  approximates  the  number  of 
cyc.es  in  the  wavelet,  vis  the  phase  shift  and  ~  is  time  shifted  by  half  the 
wavelet  auratior,. 


A  unit  amplitude  wavelet  with  80  Hz  predominant  frequency,  zero  phase 
shift  and  a  =  4  is  used  to  define  the  input  wave  fields.  A  reference  model 
with  the  same  wavelet  input  to  a  halfspace  (199  m  by  199  m  with  no  cutouts) 
provides  the  incident  field.  The  diffracted  field  is  obtained  by  subtracting 
the  incident  from  the  full  field.  Output  consists  of  time  histories  of  the 
diffracted  field  at  points  around  the  vertex  and  snapshots  of  nodal  vel¬ 
ocities  at  selected  times. 

2.3  DIFFRACTION  FROM  PLANE  SH-WAVES 

SH-waves  are  relatively  unimportant  in  mesa  applications  compared  to 
P-SV-waves  because  surface  sources  are  seldom  of  SH  type.  Yet  scalar  SH- 
waves  have  the  advantage  that  analytic  diffraction  solutions  are  available, 
whereas  no  solution  is  available  for  P-SV  diffraction.  To  investigate  val¬ 
idity  of  the  finite  element  solutions,  especially  effects  of  the  stepped  in¬ 
terface  produced  by  square  elements,  diffractions  from  SH-waves  are  studied 
and  compared  with  the  analytic  solution. 

The  closed-form  solution  for  SH  incidence  to  wedges  i:  available  from 
Keller  and  Blank  (1950)  in  terms  of  algebraic  and  trigonometric  functions 
of  radial  distance  from  the  vertex,  wave  speed  and  time.  Keller  and  Blank 
use  a  step  function  as  the  input  wave.  To  obtain  the  diffracted  field  for 
a  wavelet  input,  the  analytic  solution  is  convolved  with  the  wavelet  time 
function. 

Figure  2-2  shows  sketches  of  wavefronts  in  the  diffracted  field  for 
several  wedge  angles.  The  calculated  field  consists  of  two  components, 
reflections  off  the  two  interfaces  and  a  circular  diffracted  wave  from  the  ver¬ 
tex.  For  certain  regions  of  the  mesa,  such  as  points  A  and  B  in  Fig.  2-2, 
the  diffracted  field  has  both  components  with  the  circular  portion  lagging 
behind.  Thus,  the  time  of  arrival  is  controlled  by  the  reflected  component, 
not  by  the  circular  component  as  intuition  would  suggest.  This  is  confirmed 
by  the  finite  element  solution. 

To  simulate  2-D  plane  SH  waves,  symmetric  boundary  conditions  are  en¬ 
forced  on  both  the  left  and  right  side.  The  diffracted  field  from  the  finite 
element  analysis  is  cross-plotted  against  the  Keller-Blank  solutions  in 
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Figs.  2-3  to  2-8  for  various  reentrant  angles.  It  is  obvious  that  the  finite 
element  solution  agrees  very  well  with  the  analytic  solution.  When  the  inter¬ 
face  is  inclined  at  angles  which  give  a  smooth  stepped  interface  (such  as  0, 

45  and  90  degrees),  the  finite  element  model  gives  very  good  results.  For 
other  angles  with  a  less  uniform  stepped  interface,  the  finite  element  model 
is  able  to  capture  the  predominant  wavelet  diffraction  but  generates  some  os¬ 
cillation  at  late  time.  The  amplitude  of  this  oscillation  appears  to  be  a 
direct  function  of  distance  from  the  stepped  interface. 

Figure  2-9  shows  snapshots  of  the  finite  element  solution  for  the  geom¬ 
etry  shown  in  Fig.  2-8.  Both  the  full  wave  field  and  the  diffracted  field 
are  plotted. 

To  investigate  the  nature  of  the  late-time  oscillation,  a  mumerical  study 
with  the  same  finite  element  model  but  a  lower  frequency  (40  Hz)  wavelet  is 
made.  The  geometry  giving  the  worst  oscillation  (Fig.  2-8)  is  used.  Figure 
2-10  shows  the  diffracted  field  for  the  lower  frequency  wavelet  input.  The 
late-time  oscillation  is  reduced  significantly. 

To  further  illustrate  the  dependence  of  the  oscillation  on  the  frequency 
content  of  the  input  wavelet,  another  computation  is  made  with  a  coarse  grid 
of  2  m  square  elements  and  40  Hz  wavelet  input.  The  frequency  resolution  of 
this  coarse  grid  model  is  about  80  Hz.  Again,  the  input  wavelet  has  a  fre¬ 
quency  of  about  half  the  resolution  of  the  grid.  The  late-time  oscillation 
is  very  evident,  as  illustrated  in  Fig.  2-11. 

Finally,  a  finite  element  model  with  the  coarse  grid  and  skewed  elements 
is  used.  The  skewed  elements  eliminate  the  stepped  interface  but  impose  severe 
time-step  size  constraints  for  interfaces  inclined  at  angles  close  to  45  de¬ 
grees.  Figure  2-12  compares  the  diffracted  field  from  the  skewed  elements  to 
the  Keller-Blank  solution.  Again,  the  late-time  oscillation  is  reduced  dras¬ 
tically.  This  clearly  shows  the  relationship  of  the  oscillation  to  the  stepped 
interface. 

This  study  shows  that  finite  element  modeling  with  square  elements  can 
capture  the  predominant  response  of  the  diffracted  field.  The  stepped  inter¬ 
face  does,  however,  introduce  some  oscillation  to  the  late-time  solution. 
Whereas  the  finite  element  grid  can  be  used  to  resolve  up  to  160  Hz  for  reg¬ 
ular  wave  propagation  analysis,  the  singularity  of  the  diffracted  field  causes 


8 


a  slight  degradation  of  the  late-time  solution  for  an  80  Hz  input.  The 
finite  element  solution  gives  better  results  as  the  frequency  of  the  input  is 
lowered.  A  lower  frequency  wavelet  has  a  longer  wave  length  and  the  incoming, 
reflected  and  diffracted  wave  trains  will  interfere  with  one  another  on  the 
time  history  plots,  making  interpretation  difficult.  All  subsequent  body  wave 
calculations  will  be  made  with  the  baseline  model  of  1  m  square  elements  sub¬ 
jected  to  an  80  Hz  wavelet  input. 

2.4  DIFFRACTIONS  FROM  PLANE  P-SV  WAVES 

The  diffracted  field  from  a  vector  (P-SV)  field  is  more  complex  than 
that  from  a  scalar  field.  Incident  field  polarized  as  pure  P-wave  (where 
the  particle  motions  are  in  the  direction  of  wave  propagation)  or  pure  SV- 
wave  (particle  motions  perpendicular  to  direction  of  wave  propagation)  gen¬ 
erates  both  P  and  S V  components  in  the  diffracted  field.  The  analysis  will 
consider  first  diffractions  for  P-wave  incidence. 

A  two-dimensional  plane  strain  model  is  used  to  simulate  planar  P-waves 
with  symmetry  conditions  imposed  on  both  the  left  and  right  boundaries.  Dif¬ 
fracted  wavefronts  are  shown  in  Fig.  2-13.  Note  that  the  diffracted  field  is 
more  complicated  than  for  the  SH  case  shown  in  Fig.  2-2.  Both  P-  and  SV-waves 
are  reflected  off  inclined  interfaces  according  to  Snell's  law.  Diffraction 
from  the  vertex  has  both  P  and  SV  circular  fronts  propagating  at  their  res¬ 
pective  wave  speeds.  Of  course,  only  the  P-wave  is  reflected  at  normal  in¬ 
cidence  (horizontal  interface).  Figures  2-14  to  2-19  are  cross  plots  of  the 
x  and  y  components  of  the  diffracted  solution.  Figures  2-20  to  2-24  show 
vector  plots  of  the  diffracted  field  for  different  geometries.  Note  that 
most  of  the  vector  plots  show  only  the  center  100  by  100  m  portion  of  the  grid. 

For  the  SV  incident  case,  asymmetric  boundary  conditions  are  used  in  the 
left  and  right  sides.  Figure  2-25  shows  wavefronts  of  the  diffracted  field. 
Only  an  SV  wave  is  reflected  for  SV  incidence  beyond  the  critical  value.  This 
is  illustrated  in  the  bottom  sketch  in  Fig.  2-25.  Dirfraction  from  the  vertex 
has  both  P  and  SV  components,  figures  2-26  to  2-31  are  cross  plots  of  x  and  y 
components  or  the  diffracted  solution.  Note  the  different  time  scale  to 
accommodate  the  slower  incident  SV-wave  and  the  new  legend  for  x  and  y  com¬ 
ponents  c *  velocity.  rigures  2-32  to  2-35  show  vector  plots  of  the  diffracted 
fipld  for  d! rrer 'r‘  ceo~etries. 


2.5  DIFFRACTIONS  FROM  RAYLEIGH  WAVES 


Rayleigh  waves  are  generated  by  applying  a  vertical  velocity  over  five 
nodes  along  the  x-axis,  see  Fig.  2-1.  The  same  two-dimensional  plane  strain 
model  is  used.  The  boundary  conditions  are  symmetry  on  the  right  and  absorb¬ 
ing  on  the  other  three  sides.  In  addition  to  Rayleigh  waves,  the  normal  vel¬ 
ocity  input  generates  P- and  SV-waves ,  although  at  a  much  lower  amplitude. 

The  P-wave  speed  of  2600  m/sec  is  much  higher  than  that  of  the  Rayleigh  wave 
(at  1440  m/sec,  96  percent  of  SV  wave  speed).  Thus,  arrival  of  P  and  Ray¬ 
leigh  waves  at  the  vertex  will  be  well  separated.  The  Rayleigh  wave  can  be 
distinguished  from  the  SV-wave  by  its  retrograde-elliptical  particle  motion. 

On  crossplots  of  x  and  y  velocity  components,  this  will  show  up  as  a  90-degree 
phase  difference  between  the  two.  In  order  to  separate  the  three  different 
wave  types  even  more,  a  wavelet  frequency  120  Hz  is  used.  The  response  in 
the  second  quadrant  contains  only  the  diffracted  solution  since  this  is  a 
shadow-zone.  Thus,  the  reference  field  is  a  200-by- 100-node  (199  m  by  99  m) 
grid  corresponding  to  the  bottom  half  of  the  baseline  model.  The  diffracted 
field  is  obtained  by  subtracting  the  reference  from  the  full  field. 

Figures  2-37  to  2-40  are  cross  plots  of  x  and  y  velocities  of  the  dif¬ 
fracted  solution.  Note  the  different  scale  used  for  this  series  of  plots. 

The  low  amplitude  early  arrivals  are  diffractions  from  the  P-wave.  The  higher 
amplitude  response  is  the  diffracted  Rayleigh  wave  with  the  90-degree  phase 
difference  between  the  two  components.  The  diffracted  Rayleigh  wave  travels 
along  both  the  inclined  free  surface  (transmitted)  and  the  positive  x-axis 
(reflected),  the  former  always  having  a  higher  amplitude  than  the  latter. 
Figures  2-41  to  2-45  show  vector  plots  of  the  diffracted  field  for  different 
geometries. 


Figure  2-3.  Comoarison  of  Keller-81ank  and  FLEX  soluti 
incident  from  below  to  90-degree  mesa. 


Figure  2-4.  Comparison  of  Keller-Blank  and  FLEX  solution  for  SH 
wave  incident  to  70-degree  mesa. 
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Figure  2-10.  Comparison  of  Keller-Blank  and  FLEX  solutions  for  geometry 
shown  in  Fig.  2-8,  40  Hz  SH  wavelet  input 


Figure  2-12.  Comparison  of  Keller-Blank  and  FLEX  solutions  for  the  geometry 
in  Fig.  2-4,  40  Hz  SH  wavelet  input,  coarse  grid,  skewed 
elements. 
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Figure  2-30.  Time  histories  of  diffracted  field  from  SV-wave  incidence  to 
90- degree  mesa,  offset  45  degrees. 
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(a)  Diffracted  field  for  center  grid  at  time  0.09238  s 
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igure  2-45.  Diffractions  from  Rayleigh  wave  incidence  to  20-degree  mesa 
showing  full  wave  field  for  the  center  position  of  the  grid 


SECTION  3 

SELF-SIMILAR  THEORY  OF  VECTOR  WAVE  DIFFRACTION 

3.1  INTRODUCTION 

This  section  develops  a  new  mathematical  theory  of  vector  wave  dif¬ 
fraction.  The  theory  is  based  on  self-similar  wave  solutions  in  the  so- 
called  generalized  diffractor.  Fig.  3-1.  This  idealized  diffractor,  des¬ 
ignated  n,  is  made  of  M  contiguous  sectors  covering  the  plane  from  a  common 
vertex.  It  simulates  a  wide  variety  of  physical  structures  diffracting 
elastic,  acoustic  or  electromagnetic  waves. 

Although  many  efforts  have  been  made  over  the  years  to  solve  vector 
diffraction  problems  in  applied  physics,  geophysics,  and  mechanics,  they 
have  all  failed.  These  failures,  however,  are  due  to  mathematical  limitations 
rather  than  physics.  In  all  cases,  the  underlying  culprit  is  essential 
coupling  of  two  or  more  wave  functions  at  the  vertex--excluding,  in  particular 
separable  solutions  of  the  wave  equations.  Therefore,  the  following  theory 
is  based  on  a  completely  general  vector  formulation  of  wave  equations  and 
boundary  conditions,  with  no  further  reference  to  physics.  Since  no  charac¬ 
teristic  length  exists  in  the  generalized  diffractor,  homogeneous  solutions 
of  degree  zero  (i.e.,  self-similar  solutions)  provide  the  basis  for  analysis. 

The  theory  is  presented  in  the  following  nine  subsections,  covering: 

1)  the  governing  vector  equations;  2)  homogeneous  vector  solutions;  3)  the 
resulting  vector  boundary  value  problem  (VBVP);  4)  characteristic  mapping 
of  the  boundaries;  5)  normal  forms  of  the  VBVP;  6)  solutions  of  all  normal 
forms;  7)  consistency  conditions;  8)  reduction  to  a  vector  integral  equation; 
and  9)  the  canonical  problem  (two  wave  functions). 

3.2  EQUATIONS  GOVERNING  THE  GENERALIZED  DIFFRACTOR 

In  two  space  dimensions,  the  generalized  diffractor,  n  consists  of  M  semi 
infinite  wedges  formed  by  M  directed  lines  emanating  from  the  common  vertex. 
These  directed  lines  (interfaces)  partitioning  the  plane  are  identified  by 
polar  angles,  L  ,  m  =  1,  M,  and  their  union  is  designated  L.  Geometry  and 
polar  coordinates  are  illustrated  in  Fig.  3-1. 

A  total  of  N  wave  functions,  f^  to  f^  are  supported  by  the  diffractor 
in  domains,  ft  ,  n  =  1,  N.  The  number  of  wave  functions  per  wedge  is  denoted 
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j  ,  m  =  1,  M,  and  by  definition  each  function  satisfies  a  scalar  wave  equation. 

m  j 

Consequently,  the  N-vector  unknown,  f  =  (f  . ..f^)  is  governed  by  the  normal 
form  of  the  vector  wave  equation  (vector  d'Alembertian) , 

B  f  =  c2Af  -  f  t  =  0  in  ft  (3-1) 

where  Af  =  f  +  f  /r  +  fQQ/r2  is  the  scalar  Laplacian  acting  on  the  elements 
rr  r  00 

of  f,  and  c  is  a  diagonal  N  x  N  matrix  of  real  constants.  Partial  derivatives 
are  of  course  indicated  by  r,  0  and  t  subscripts.  This  diagonal  form  of  the 
operator,  B  does  not  preclude  off-diagonal  terms  in  c,  only  requiring  that  the 
fuller  matrix,  e.q.  c,  be  diagonalized  by  a  linear  transformation,  S  as  c  = 
S-1cS,  whence  the  solution  is  f  =  Sf. 

The  N  functions  defined  in  M  wedges  are  coupled  by  boundary  conditions  on 
set  the  union  of  Tl-domain  limit  points.  In  addition  to  interfaces,  L, 
these  points  include  the  limiting  contour,  C  =  Cq  +  Cm  of  circles  bounding 
the  vertex  and  infinity.  The  union  of  interfaces  and  contours  is  illustra¬ 
ted  in  Fig.  3-2.  General  boundary  conditions  on  the  limit  points  are  written 


rJr  +  Tft  +  nf0/r  =  g  on  an  =  L  +  C  (3-2) 

where  c,  x  and  n  are  linearly  independent  (disjoint)  N  x  N  matrices  and  g  is 
an  N-vector. 

The  matrix  coefficients  and  inhomogeneous  vector  in  (3-2)  are  actually 
distributions  on  3 n  and  represent  sums  over  the  union  of  contours,  namely 


M 

C  =  Z 
m=l 


m 


M 
=  Z 
m=l 


m 


M 

=  I  n, 
m=l 


m 


M 

=  Z  g, 

m=l 


m 


Terms  in  each  sum  are  linearly  independent  by  virtue  of  their  definition  on 

disjoint  boundaries,  with  the  common  vertex  point  excluded  by  contour  Cg. 

On  interface  L  ,  condition  (3-2)  prescribes  j  +  j  ,  equations  in  the 
m  r  ^m  m-1 

same  number  of  wave  functions.  On  arc  C  ,  (3-2)  prescribes  jm  equations  in 

m  m 

j  unknowns.  The  nonzero  row  and/or  column  elements  in  t  ,  t  ,  n  and  g 
vrn  m  nr  'm  m 

correspond  to  positions  o^  the  coupled  functions  in  N-vector,  f. 


3.3  SELF-SIMILARITY  AND  HOMOGENEOUS  VECTOR  SOLUTIONS 


Because  the  di^ractor  lacks  a  characteristic  length,  self-similar  solu¬ 
tions  of  the  govern' ng  equations,  (3-1)  and  (3-2),  are  admissible.  The  most 
tractable  and  directly  useful  of  these  are  homogeneous  solutions  of  degree 


zerr. 


;t  lcn:.  where  '(r.h,!)  ~  t'" r(r/ 1 ,  ■> ) . 


!n  or ’ nc ’ o 1 c ,  homogeneous 
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solutions  of  non-zero  degree  can  be  obtained  by  differentiation  or  integra¬ 
tion  of  the  degree  zero  solution. 

Introducing  the  scalar  similarity  variable,  p=  r/t  into  the  wave  equation, 
(3-1)  using  the  relation,  rf^  =  -tft  =  pf  ,  gives 

p2(c2-p2I)f  +  c2faa  +  p(c2-2p2I )f  =  0  (3-3) 

pp  00  p 

where  I  is  the  N  x  N  identity  matrix  (ones  on  the  main  diagonal  and  zeros  off). 

The  most  effective  approach  to  solving  this  self-similar  form  of  the  vector 
wave  equation  is  to  transform  it  to  characteristic  coordinates.  The  charac¬ 
teristic  equation  is  readily  found  to  be 

p2(c2-p2I) (de)2  +  c2(dp)2  =  0 

and  solving  for  d0/dp  gives  the  differential  equation 

-  f  I  -  6t  •  -  iifejj'Vp) 


The  characteristic  function,  0+  has  a  pole  at  p  =  0,  branch  point  at  pi  =  c, 
and  a  double  zero  as  p-*».  Principal  values  of  the  square  root  and  inverse 
cosine  functions  are  assumed,  while  multivaluedness  is  included  explicitly  with 
the  ±  sign.  Multiplying  by  dp  and  integrating  yields 

z+  =  01  ±  cos”lc/p  =  x  ±  y  (3-4) 


where  the  new  coordinate  matrices,  z+  =  x  ±  y  are  constants  of  integration. 

Transforming  the  similarity  form,  (3-3)  to  characteristic  coordinates 
by  replacing  the  derivatives  as 

_  ’  fp  =  Mz 

yields  the  normal  vector  form. 

f  =  0 


v  vv 


+  B-fz 


(3-5) 


z  +  z- 
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which  of  course  can  be  integrated  directly  giving  the  classic  solution 

f(p,e)  =  F+(z+)  +  F_  (z_  ) .  (3-7) 


3.4  THE  VECTOR  BOUNDARY  VALUE  PROBLEM 

Particular  solutions  for  F+  are  determined  from  the  boundary  condition, 


(3-2).  Introducing  the  similarity  variable,  p  =  r/t  gives 

nfg/r  +  (e-px)f  /t  =  g 

Therefore,  to  admit  homogeneous  solutions,  matrices  c,  n  and  t  must  have 
the  same  degree  of  homogeneity,  h  say,  while  vector  g  is  necessarily  of 
degree  h-1.  This  means  that  g(r,t)  =  t ^  ^g(p)  and  e.g.,  n(r,t)  =  t  n(p)> 
whence  dividing  the  above  by  t^1”^  yields  the  boundary  condition, 

n/pf0  +  (?-PT)f  *  g(p) 

Substituting  (3-7)  into  (3-5),  the  partial  derivatives  of  f  become 

fe  =  f;  +  K  ■  V  e/V  +  8-FJ 

Since  3  =  -3+  these  can  be  written 

F; +  F: =  fe 
f;  -  f;  =  3;xfp 

Therefore,  defining  densities,  d+  such  that  on  the  boundary 

d+  ~  f6  ’  d-  E  B+lfp 

yields  the  following  boundary  value  problem  on  3J1 

f;  ±  k  » d± 

n/pd  +  +  ( C-  P"t )  B  +c’_  =  g 


(3-8) 


(3-9) 


(3-10) 


(3-11) 

(3-12) 


for  determination  of  F!  in  ft.  Because  d,  and  d  (i.e.,  f.  and  f  )  are  not 

±  '  “  Dp 

independent,  the  multivalued  problem,  (3-11)  represents  a  consistency  condition 


between  f.  and  f  . 

e  p 


3.5  CHARACTERISTIC  MAPPINGS 

The  analytic  functions,  F+(z+)  are  determined  by  solving  the  boundary 
value  problems  defined  in  (3-11)  and  (3-12)  on  the  limit  set,  311.  The 
characteristics  map  31!  to  3ft,  a  semi- infinite  rectangle  in  the  z+(or  z_) 
plane  enclosing  the  wave  domain,  ft.  Recall,  because  z+  is  a  diagonal  N  x  N 
matrix  of  characteristics,  this  mapped  domain  actually  represents  N  complex 
planes. 

The  z+  plane  is  illustrated  in  Fig.  3-3,  and  z_  is  a  simple  reflection 
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about  the  x  (real)  axis.  The  bounds  on  x  are  «><x<  t>+  u>  ,  where  d  and  w  are 

diagonal  N  x  N  matrices.  Elements  of  d  define  angles  of  the  lower  interfaces 

and  those  of  u  define  wedge  angles,  in  the  N  wave  domains.  In  particular, 

wedge  m  supports  j  wave  functions,  then  for  each  of  the  jm  functions 

j  e  L  and  u>  =  l  -  L  .  The  bounds  on  y  are  -<=°<y<Tr/2,  where  characteris- 
n  m  n  m+1  m  ~  - 

tics  are  real  for  y>0  (hyperbolic  subdomain)  and  complex  for  y<0  (elliptic 
subdomain).  The  circular  arc  at  infinity,  maps  to  the  end  of  the  rectangle 
at  y  =  u/2 ,  while  the  arc  around  the  vertex,  CQ  goes  to  the  opposite  end, 
at  y  -*•-  ». 

The  boundary  value  problem  for  F+  is  complicated  by  the  corners  on 
and  mixed  behavior  of  the  functions,  i.e.,  hyperbolic- elliptic.  These 
difficulties  can  be  eliminated  by  a  suitable  transformation  of  the  characters 
tics.  Consider  a  succession  of  mappings  beginning  with  the  linear  trans¬ 
formation, 

ITU)  1(Z+-J?) 

which  scales  the  real  coordinate  of  the  wave  domain  between  0  and  tt.  Taking 
the  cosine, 

cos(iruf  x(z+-  «?)) 

maps  the  boundary  to  the  real  axis,  with  z+  and  z_  domains  going  to  the  lower 
and  upper  halfplanes,  respectively.  The  hyperbolic  subdomain  closure  maps 
along  characteristics  to  the  real  axis  between  -I  and  +1.  Finally,  taking 
the  reciprocal, 

w  =  [cos(iru)" 1  (z+-t?) )  ]~ 1  =  u  +  iv  (3-13 

maps  the  elliptic  z+  subdomain  to  the  upper  halfplane,  z  to  the  lower  and  the 
wedge  vertex  to  the  origin.  The  hyperbolic  subdomain  maps  entirely  to  the 
real  axis,  on  - I >u> I . 

The  conformal  mapping,  w(z+)  is  illustrated  in  Fig.  3-4  and  maps  each 

interface,  lm  to  segments  of  the  real  axis  in  jm+jn)_ ^  of  the  N  wave  domains. 

Endpoints  of  L,  at  p  =  0,°°  go  to  w  =  0,+u^,  where  u^  =  [cos(tt2u)"1/2)]‘1  . 

For  the  j  wave  functions  defined  in  wedge  m  (above  interface  m),  L  maps 
m  m 

to  [0,+u^],  while  for  the  jm_^  functions  in  wedge  m-1  Lm_^  maps  to  [0,-u^]. 
These  segments  [0,±uj  extend  from  zero  to  the  right  or  left  respectively  and 
typically  wrap  around  ±°°. 


3.6  REDUCTION  OF  F^  ±  F|  =  d+  TO  NORMAL  FORMS 

Since  densities,  d,  =  fQ  and  d  =  6*,1f  are  not  independent,  the  multi- 

+  D  -  T  P 

valued  vector  boundary  value  problem,  (3-11)  must  first  be  solved  for  their 
relationship.  Here  they  are  reduced  to  certain  normal  forms,  which  will  next 
be  solved  in  general.  By  virtue  of  the  conformal  mapping,  multivalued  sub¬ 
scripts  on  F'  now  correspond  to  limits  from  above  and  below  the  real  axis 
in  the  w  =  u  +  iv  plane,  i.e. 

F;(u)  =  lim  F'(w), 
v->0± 

rather  than  characteristics.  By  contrast,  subscripts  on  d  merely  distinguish 
righthand  sides  for  the  two  boundary  value  problems. 

The  boundary  value  expression  given  by  the  minus  sign  in  (3-11)  defines 
the  normal  form. 


The  remaining  expression, 

F1  +  F'  =  d 

can  be  reduced  to  normal  form  by  introducing  matrix  functions,  D+  which 
satisfy  the  homogeneous  problem, 

D+  +  D_  =0  (3-15) 

Multiplying  by  D+ and  replacing  D+F_)  by  -D_F_)  through  (3-15)  gives 

D+F+-  D_F;  =  D+d+  (3-16) 

Finally,  the  boundary  value  problem  for  D+,  (3-15)  is  reduced  to  normal 
form  by  rewriting  it  as 

D+DI 1  =  D'M)+  =  -I 

Because  D+  and  D~x  commute, taking  the  principal  value  logarithm  yields 

LnD  +  "  LnD_  =  i ttI  (3-17) 

Thus,  the  original  expression,  (3-11)  is  equivalent  to  the  three  normal 
forms,  (3-14),  (3-16)  and  (3-17). 
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3.7  SOLUTION  OF  THE  NORMAL  FORM,  B+  -  =  b 

The  normal  boundary  value  problem  is  defined  by  (3-14),  rewritten  in 
general  as 


B+  "  B_  =  b 


(3-18) 


To  solve  for  B(w),  the  inhomogeneous  vector  or  matrix,  b  is  replaced  by  an 
integral  over  an  using  the  sifting  property  of  the  delta  function,  i.e., 

b(u)  =  /ds<$(s-u)b(s) 

an 

Factoring  the  diagonal  matrix  of  delta  functions,  6(s-u)  as 
fi(s-u)  =  ^rlim[s-wp  -  J^-limCs-w]-1 


gives 


where 


b(u)  =  lim  E[w;b]  -  lim  E[w;b]  =  E ,  -  E 


£[w;b]  =  2~j  /ds[s-w]-1b(s) 

'an 


The  boundary  values  of  E  are  found  to  be 


E±[u;b]  =  b(u)  +  Ep(u;b) 


(3-19) 


(3-20) 


where  Ep  is  the  principal  value  integral ,  i  .e. ,  E[u;b]  in  (3-19)  with  the  point 
s  =  u  excluded.  Clearly  a  particular  solution  for  B(w)  is 

B(w)  =  Z[w;b]  (3-2i; 


3.8  THE  CONSISTENCY  CONDITION  ON  DENSITIES 

Normal  forms  of  the  multivalued  boundary  value  problems  on  F|  can  now 
be  solved.  From  (3-14)  the  sc’ution  in  terms  of  d  is 

F;  -  F;  =  d_  =»  F'(w)  -  E[w;d_]  (3-22; 


while  for  d+,  (3-16)  and  (3-17)  yield 


f;  +  f;  =  d+ 


F'(w)  =  D"l(w)E[w;D  d^l 


(3-23) 


am  a  n*  r*  v  jr*  w~±  w~*  r*  »y  ^ 


J  .  t.  ,,,  ,J  ,  ,1,1 „  t,  y  w  q  „  m,  . . .,.  ^ 


i  • » w  «■  i  ■  r»  www  nr 


Boundary  values  of  these  solutions  are 

(  ±|d  +  £_  [u ;  d  ] 

%<“>  =  M  -  ^  ' 

l  K  ±  d-^U-.d^ 


(3-24) 


which  must  be  equal,  whence,  equating  singlevalued  parts  and  multivalued  parts 
gives  the  consistency  conditions, 

|d+  =  Ep[u;dJ  (3-25) 

|d_  =  D“+lZ  [u;D+dJ  (3-26) 

This  integral  transform  pair  provides  the  relation  between  densities  d+  on 
the  boundary.  Eliminating  d+  or  d_  yields  an  integral  equation  for  each  as 

£l+  =  IptujD^ZptsjD^^]  (3-27) 

^d_  =  D^1Ep[u;D+2p[s;d_  ]]  (3-28) 

To  complete  these  expressions,  (3-17)  is  solved  for  LnD(w)  using  (3-21), 
as 

LnD(w)  =  E[w,iTrI]  =  j  /ds[s-w]_1 

Taking  the  matrix  exponential  gives 

D(w)  =  eilrZ^3  (3-29) 

where  l[w]  =  E[w;  I],  and  from  (3-20),  boundary  values  are 

0>)  -  eiw(±I/2  +  2PM)  .  ±(eii>£p[u3 

3.9  REDUCTION  OF  THE  BOUNDARY  VALUE  PROBLEM  TO  AN  INTEGRAL  EQUATION 

To  solve  the  boundary  condition  equation,  (3-12),  d+  or  d  is  replaced 
through  (3-25)  or  (3-26).  Since  the  normal  form  of  solution  for  F'(w) 
is  given  by  (3-22)  in  terms  of  d_ ,  it  is  natural  to  eliminate  d+.  Therefore, 
inserting  (3-25)  in  (3-12)  gives  the  integral  equation. 

2n/p£p[u;d_]  +  (c-px)B+d_  =  g  (3-31) 


63 


However,  in  order  for  this  to  be  self-consistent,  d_  must  also  satisfy  the 
consistency  integral  equation,  (3-28).  These  equations  are  combined  by 
operating  on  (3-31)  with  EpCujpD.p'"1] ,  whence 

2£p[u;D+Ep[u;d_]]  +  Ep[u;  pD+rf3(c-pT)B+d_]=  2p[u;pD+rf  lg] 

and  replacing  the  double  integral  by  (3-28),  giving 

^0+d_  +  Ip[u;pD+rf  Mc-p-Oe+dJ  =  Zp[u;  pD+n"  xg3  (3-32) 

This  equation  is  valid  provided  rf  1  exists.  If  n  is  not  invertible,  the 
implication  is  that  elements  of  vector  d__  are  linearly  dependent,  i.e., 
some  degree  of  degeneracy  exists.  To  proceed  with  this  singular  case,  however, 
physics  must  be  invoked  to  evaluate  admissible  forms  of  n»  t,  and  g. 

The  single  integral  equation,  (3-32)  is  the  final,  self-consistent 
expression  of  the  boundary  value  problem,  (3-11)  and  (3-12).  By  virtue  of 
its  N-vector  unknown  and  dependence  on  N  complex  variables,  no  standard 
analytical  reduction  of  the  problem  appears  possible.  Thus,  to  determine 
F'(w)  it  is  necessary  to  solve  (3-32)  for  d_  on  3ft  and  evaluate  E[w;d_]. 

In  full,  this  gives  F'(w)  as 

F'(w)  =  E  [u;d_]  =  ptr  /ds[s-w]_1d  (s)  (3-33) 

p  CTU  ^ 

by  solving 

jD+d_  +  Ep[u;v-1pd_]  =  Ep[u;v_1g]  (3-34) 

where 

V  =  nD"+l/p  ,  D+  =  ieilTEP[u] 

y  =  (c-px)B+  ,  B+  =  c(p/pTFc7)  (3-35) 

Direct  integration  yields  F(w)  as 

F<w)  =  3^dsLr'l:5-w3d.{s)  (3.36) 

where  Ln[s-w]  is  a  diagonal  matrix  of  logarithms  and  initial  conditions  are 
specified  on  the  mapping  of  C  to  3ft. 


3.10  THE  CANONICAL  PROBLEM,  N=2 

The  generality  of  this  formal  solution  for  F'(w),  (3-33)  in  terms  of 
a  singular  integral  equation  precludes  further  resolution  of  the  actual 
computational  problem.  To  proceed,  it  is  necessary  to  consider  a  specific 
class  of  problem  solutions.  In  particular,  specializing  to  the  lowest  order 
vector  problem,  N=2  yields  the  (binary)  canonical  diffraction  problem.  This 
involves  only  two  wave  functions,  either  with  both  supported  in  one  wedge, 

M=l,  or  one  supported  in  each  of  two  wedges,  M=2.  Thus,  the  binary  problem 
is  in  terms  of  2-vectors  and  2x2  matrices  only. 

Rather  than  consider  specific  boundary  conditions,  i.e.,  values  of 
n>C  and  x  in  (3-2),  because  they  are  disjoint  and  2x2,  their  general 
structure  can  be  determined.  In  particular,  there  are  only  three  possible 
pairs  of  disjoint  2x2  matrices,  namely 

(ID  ,  (?J) 

(ID  .  (!!) 

(ID  ,  (SI) 

The  first  pair  is  invertible  (nonsingular),  while  the  other  two  are  not 
(singular).  Consideration  of  typical  boundary  conditions  in  physical 
problems  show  that  the  first  pair  of  matrices  corresponds  to  the  M=1  case, 
the  second  pair  to  M=2,  while  the  third  corresponds  to  a  degenerate  scalar 
case  which  will  be  neglected.  The  canonical  coefficients  can  therefore 
be  written, 

M=l;  v  =  a(p)^J  °)  ’  y  =  b(p)(i  j)  (3-37) 

M=2;  v  =  a(p)  ,  y  =  ^  ^  b(p)  (3-38) 

where  a(p)  and  b (p)  are  2x2  diagonal  matrices  which  depend  on  the  specific 
application. 
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The  M=1  case  satisfies  the  assumptions  made  in  deriving  the  integral 
equation,  (3-34),  hence,  substituting  (3-37)  into  (3-34)  yields  the  2-vector 


equation 


.  E  [uA/aiV 

(b2/a2)dl/J  p  \92/a2/' 


(3-39) 


where  subscripts  on  D  and  d  are  dropped  for  convenience.  The  position  reversal 
of  d1  and  d^  in  the  integral  operator  precludes  further  simplification  (with¬ 
out  introducing  a  multiple-integral  scalar  equation).  However,  the  final 
system  of  equations  is  solvable  using  successive  approximations.  Writing 
it  as 


2D' 


'oya^d. 

PLU  *V (b27 a2 ) d 


(3-40) 


and  assuming  starting  values  for  d^  and  dg,  this  recurrence  relation  yields 
a  new  approximation  for  d^  and  d^^  respectively.  Repeating  the  process  yields 
successively  better  approximations  for  vector  d. 

The  singular  M=2  case  does  not  satisfy  the  assumption  leading  to  the 
formal  solution  (3-34),  i.e.,  existence  of  v"1.  It  is  necessary  therefore 
to  directly  examine  boundary  condition,  (3-12),  rewritten  as 


vd+  +  yd_  =  g. 

Substituting  (3-38)  gives  two  equations  relating  components  of  d+  and  d_, 

aldl_  +  a2d2-  *  91 

bldl+  +  b2d2+  s  S2  l3‘‘ 


Eliminating  d+  in  (3-41)  through  (3-25)  gives  the  scalar  equation 

bl*Vul’dl-J  ”  b2^p^u2’al^a2dl-^  ~  9  2^  ~  b2^'p^u2’^l'/a2^  (3-42) 


while  eliminating  d_  through  (3-26)  gives 

aj/DjEp[ Uj ,D|di+)  -  a 3/ ( u2 » ^2 ( ^1^ bg ^ d ^+]  =  g^/2  -  a  2^  2^  2^  2) 

These  integral  equations  can  be  solved  for  d^_  and  d^+  respectively,  and  the 
other  components  follow  from  (3-41).  Since  the  resulting  solutions  for  d+ 
and  d_  are  not  independent, clearly  one  or  the  other  is  redundant.  This 
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can  only  be  the  case  if  components  of  the  g  vector  are  depenc 
c)2  =  92^l)  *  therefore,  in  the  M=2  case  (for  which  v"1  and  y 
this  singular  problem  can  only  be  solved  if  g  is  a  degenerate 
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Figure  3-3.  Mapping  of  the  wave  domain,  ft  into  the  z 
by  the  characteristic  transformation. 


W  PLANE 


Mapping  of  the  wave  domain, ft  to  the  upper  half 
plane  by  conformal  transformation,  w(z,). 


Figure  3-4 


SECTION  4 


DISCUSSION  AND  CONCLUSION 


4.1  OVERVIEW 

This  report  has  addressed  numerical  and  analytical  aspects  of  wave 
diffraction  from  corners  common  to  geologic  and  topographic  structures.  The 
numerical  models  examined  are  limited  to  reentrant  corners  (angles  greater 
than  180  degrees);  however,  the  mathematical  analysis  is  completely  general. 

The  numerical  part  has  examined  accuracy  and  sensitivity  of  numerical 
diffraction  to  di scretization  errors.  However,  a  quantitative  study  of 
accuracy  was  only  feasible  for  the  relatively  unimportant  scalar  case  of  SH- 
wave  diffraction.  Therefore,  as  a  first  major  step  in  evaluating  accuracy 
of  P-SV  numerical  diffractions,  the  underlying  mathematics  were  studied. 

The  resulting  analytical  part  of  the  report  provides  a  new  mathematical 
basis  for  wave  diffraction  from  a  generalized  model. 


4.2  NUMERICAL  DIFFRACTION  EXPERIMENTS 

It  is  necessary  to  study  elastic  corner  diffraction  using  numerical  meth¬ 
ods,  since  no  analytical  solutions  are  available.  However,  there  are  questions 
regarding  accuracy  which  must  first  be  answered.  Discrete  numerical  wave 
solutions,  using  finite  differences  or  finite  elements,  implicitly  assume  regu¬ 
larity  of  the  solution.  When  initial  data,  boundary  data  or  boundary  geom¬ 
etry  are  discontinuous,  the  exact  solutions  may  exhibit  singularities  which 
are  excluded,  i.e.,  not  represented  adequately,  in  the  numerical  solution-- 
without  special  fixes.  This  is  particularly  true  at  corners,  where  the  tangent 
exhibits  a  jump  rather  than  turning  continuously,  especially  when  the  corner 
is  reentrant.  This  problem  can  be  clarified  by  comparing  finite  element  dif¬ 
fractions  to  existing  scalar  diffraction  solutions. 

The  results  of  Section  2.3,  comparisons  of  the  Keller-Blank  scalar  dif¬ 
fraction  solutions  to  SH  finite  element  solutions,  show  that  simple  finite 
element  corner  models  give  good  results  for  the  predominant  diffracted  signal. 
However,  the  good  correlation  deteriorates  slightly  at  later  time  due  to  dis¬ 
cretization  errors  associated  with  model  roughness  on  the  sloping  wedge  face. 
These  errors,  although  insignificant  for  this  early  time  study,  can  be  further 
reduced  by  more  refined  gridding  or  skewed  elements,  as  demonstrated.  There¬ 
fore,  these  experiments  indicate  that  with  no  special  f’xes,  finite  elements 
give  excellent  results  for  the  simolest  case  of  scalar  wave  diffraction. 
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The  more  realistic  case  of  P-SV  diffraction  is  much  more  complicated 
than  the  scalar  case  since  two  particle  velocity  components  are  necessary 
to  describe  the  field  at  each  point.  Also,  the  reentrant  corner  produces 
a  stress  singularity  at  later  time,  and  excites  Rayleigh  surface  waves  in 
addition  to  diffracted  body  waves.  Since  there  are  no  corner  solutions  for 
quantitative  comparison,  finite  element  simulations  of  P-SV  wave  diffraction 
were  simply  examined  for  consistency.  One  troublesome  aspect  was  graphical 
display  of  the  vector  field.  Time  histories  and  surface  plots  of  each 
component  were  used  but  the  most  successful  display  approach  was  the  simple 
vector  plot.  Vector  plots  of  the  velocity  field  will  be  used  extensively 
in  future  two-dimensional  simulations. 

The  final  group  of  calculations  examined  Rayleigh  surface  wave  diffraction 
from  the  corner.  This  is  a  particularly  interesting  physical  problem  and  has 
been  the  subject  of  numerous  experiments  and  approximate  analyses.  Rayleigh 
wave  diffraction  is  even  more  complicated  then  P-SV  diffraction  because  of 
precursor  body  waves,  as  well  as  the  surface  wave's  underlying  structure. 
However,  the  finite  element  calculations  are  no  more  complicated  than  for  plane 
P-  or  SV-wave  incidence--actual ly  simpler  because  the  Rayleigh  source  is  a 
pressure  pulse  acting  on  only  part  of  the  free-surface.  Among  other  results, 
these  calculations  show  that  surface  wave  transmission  is  always  much  stronger 
than  reflection  back  along  the  free-surface,  even  for  a  vertical  face. 
Furthermore,  the  SV  body  wave  diffracted  into  the  mesa  shadow  zone  is  quite 
strong  for  reentrant  corners. 

The  complete  set  of  finite  element  calculations  were  performed  on  a 
minicomputer  (Prime  750).  Model  size  was  typically  less  than  50,000  elements 
due  to  limited  fast  memory  size.  The  number  of  available  elements  puts 
physical  limits  on  frequency  resolution  of  the  model.  This  was  seen  in  the 
resulting  field  plots  as  fairly  long  wavelet  lengths,  compared  to  model  size. 
Some  details  of  interest  were  therefore  missed  due  to  overlapping  of  incident 
and  reflected  or  diffracted  wavelets— in  particular,  details  near  the  Rayleigh 
wave  and  diffracted  SV  wave  at  the  surface.  This  lack  of  high  resolution 
also  tends  to  produce  poor  seismograms,  since  separate  arrivals  overlap  and 
interfere.  The  solution  to  this  problem  is  to  perform  the  calculations  on  a 
supercomputer.  However,  for  purposes  of  this  study,  the  convenience 
of  in-house  computing  and  postprocessing  outweighed  the  inconvenience  of 
accessing  a  supercomputer  on  relatively  low  data  rate  telephone  lines. 


i'  m  m  h  m* 
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4.3  THEORY  OF  VECTOR  WAVE  DIFFRACTION 


The  major  contribution  of  this  report  is  the  new  mathematical  theory  of 
diffraction  derived  in  Section  3.  Because  the  analysis  is  completely  general, 
in  terms  of  M  wedge  domains  and  N  wave  functions,  irrespective  of  physics, 
the  results  apply  to  a  wide  variety  of  unsolved  wave  diffraction  phenomena. 
These  include  optical  and  other  electrodynamic  phenomena,  as  well  as  the 
baser  mechanical  phenomena  of  elastodynamics  and  acoustics--to  which  this 
report  is  addressed. 

The  mathematical  problem  of  vector  wave  diffraction,  where  two  or 
more  wave  functions  are  coupled  at  the  generalized  diffractor,  has  been 
perplexing  and  frustrating  analysts  for  well  over  100  years.  The  fundamental 
difficulty  is  essential  coupling  of  the  wave  functions  at  the  diffractor.  As 
a  consequence,  separable  solutions  of  the  wave  equations,  the  basis  for  much 
of  classical  wave  theory,  fail  to  satisfy  the  boundary  conditions  there. 

Non-separable  forms  based  on  self-similarity,  i.e.,  homogeneous  solutions, 
are  well-known  and  have  been  applied  successfully  to  scalar  problems--but  not 
to  the  vector  case.  The  failure  is  due  to  an  incomplete  vector  formulation  of 
the  problem,  rather  than  any  underlying  inconsistency  in  tie  mathematics. 
Therefore,  it  is  necessary  to  completely  vectorize  the  mathematics  by  ex¬ 
pressing  the  governing  partial  differential  equations  as  a  matrix  differential 
operator  acting  on  a  vector  unknown,  and  logically  extending  the  concepts  of 
scalar  boundary  and  initial  conditions  to  the  vector  case. 

With  the  problem  vectorized,  solution  of  the  vector  wave  equation 
oroceeds  in  much  the  same  way  as  for  the  scalar  case.  It  starts  by  trans¬ 
formation  to  characteristic  coordinates,  which  is  equivalent  to  diagonalizing 
the  matrix  partial  differential  operator,  and  ends  by  a  trivial  integration, 
yielding  the  solution  as  the  sum  of  two  characteristic  functions,  i.e.,  func¬ 
tions  of  the  two  characteristic  matrices. 

The  real  problem  then  begins,  namely,  finding  particular  forms  of  the 
characteristic  functions  satisfying  vector  boundary  conditions  on  the  union 
of  interfaces.  These  conditions  yield  a  boundary  value  problem  involving  N 
unknown  functions  in  2N  characteristic  variables,  defined  on  limit  points  of 
the  N  wave  domains. 
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To  express  this  boundary  value  problem  in  a  form  amenable  to  analysis, 
it  is  necessary  to  further  transform  the  wave  domains  so  they  cover  the  upper 
half plane,  with  all  boundaries  mapping  to  the  real  axis.  This  is  called  the 
characteristic  mapping,  because  it  turns  out  that  real  characteristics  map 
entirely  to  the  real  axis,  effectively  vanishing,  except  for  a  residual  dis¬ 
tribution  on  the  axis.  Complex  characteristics  then  become  the  only  spatial 
variables  in  the  problem,  whence  the  2N  characteristic  variables  reduce  to 
N  complex  variables  in  the  upper  halfplane,  and  their  N  complex  conjugates 
in  the  lower  halfplane. 

Thus,  the  vector  boundary  value  problem  assumes  its  fundamental  form, 
i.e.,  given  a  vector  of  N  analytic  functions  in  N  complex  variables,  find  a 
particular  solution  satisfying  boundary  conditions  on  the  real  axis.  Part  of 
this  fundamental  form  is  expressed  as  three  so-called  normal  boundary  value 
problems,  in  order  to  enforce  consistency  between  real  and  imaginary  parts  of 
the  N-vector  analytic  function  on  the  real  axis.  The  remaining  part  of  the 
problem  is  the  vector  boundary  condition  equation. 

By  virtue  of  their  definition  on  the  real  axis,  the  normal  forms  are 
solved  essentially  by  inspection,  using  the  sifting  property  of  the  delta 
function.  This  yields  all  normal  solution  in  terms  of  a  diagonal  singular 
integral  operator  in  the  N  complex  variables.  The  three  normal  problems 
expressing  the  consistency  condition  yield  an  integral  transform  pair  re¬ 
lating  real  and  imaginary  parts.  The  pair  can  be  reduced  to  a  double  integral 
equation  by  eliminating  either  the  real  or  imaginary  part  of  the  analytic 
vector  function.  This  yields  either  part  in  terms  of  a  double  singular  inte¬ 
gral  of  itself. 

Finally,  to  solve  the  original  boundary  condition  equation,  the  real  or 
imaginary  part  is  eliminated  via  the  integral  transform  pair.  Operating  on 
the  resulting  integral  equation  with  a  suitable  integral  operator,  and  re¬ 
placing  the  double  singular  integral  yields  a  standard  singular  integral 
equation.  This  form  of  the  boundary  value  problem  serves  as  the  calcula- 
tional  basis  for  numerically  evaluating  vector  wave  diffraction  solutions. 

The  form  depends,  however,  on  invertibil i ty  of  two  matrix  coefficients  in 
the  boundary  conditions.  When  either  matrix  is  non-invertible,  i.e.,  singular. 
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the  problem  is  degenerate  and  physics  has  to  be  invoked  to  give  a  consistent 
integral  equation. 

Rather  than  end  the  analysis  at  this  level  of  generality,  it  was  use¬ 
ful  to  specialize  the  results  to  the  N  =  2  canonical  problem.  This  so- 
called  binary  problem  involves  two  wave  functions  in  a  single  wedge,  M  =  1, 
or  one  wave  function  in  each  of  two  wedges,  M  =  2.  The  M  =  1,  problem  cor¬ 
responds  to  the  fundamental  P-SV  elastodynamic  problem  examined  numerically 
in  Section  2.  The  M  =  2  problem  would  be  obtained  by  adding  an  SH  wave 
function  to  the  void  wedge  in  the  model  used  to  verify  numerical  solutions 
of  the  scalar  wave  equation  in  Section  2.  The  M  =  1  problem  is  of  principal 
interest  in  this  work,  and  is  solved  rigorously  by  the  above  mathematical 
theory  since  the  boundary  matrices  are  invertible. 


4.4  CONCLUSIONS 

A  number  of  significant  conclusions  can  be  drawn  from  this  work.  First, 
regarding  the  underlying  geophysical  problem  of  diffraction  into  shadow  zones 
by  reentrant  corners  in  highland-lowland  topography--the  strongest  diffrac¬ 
tions  from  the  corner  are  due  to  the  Rayleigh  wave  incidence.  In  addition, 
the  observed  strong  transmission  of  Rayleigh  waves  around  corners  yields 
secondary  diffractions  as  subsequent  corners  are  encountered. 

Second,  regarding  finite  element  diffraction  calculations--simple 
Cartesian  finite  element  models  appear  adequate  near  the  corner  for  early 
time,  and  also  near  the  wavefront,  but  not  for  later  time  (quasistatic) 
because  of  an  incomplete  mathematical  formulation  there,  i.e.,  no  explicit 
representation  of  stress  singularities  in  the  element  shape  function.  Ana¬ 
lytical  solutions  of  P-SV  diffraction  are  necessary  to  completely  validate 
these  results,  however. 

Third  and  most  important,  regarding  the  new  mathematical  theory  of 
diffraction--the  analysis  yields  a  rigorous,  formal  solution  of  diffraction 
by  the  generalized  diffractor  and,  in  particular,  yields  a  solution  for  body 
wave  diffraction  by  a  single  wedge.  This  provides  the  future  basis  for  val¬ 
idating  the  P-SV  numerical  diffractions  described  above,  as  well  as  solutions 
to  a  number  of  interesting,  unsolved  problems  in  mathematical  physics. 
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